
NUM_MC = 10000
set.seed(0)

# SIMULATION PARAMETERS ################################################
N = 75
T = 75
b = matrix(c(0.1,0.2,0.1,0.1,0.05,0.05),,1)
#b = matrix(c(0,0,0,0,0.25,0.25),,1)
DGP_AR = TRUE ### AR if TRUE; MA if FALSE ###
rho = .25 ### In case of AR ###
q = 10
MAcoef = matrix(c(rep(1/q,q)),,1) ### In case of MA ###
true_beta = matrix(c(.1,.1),2,1)
binom_N = 5
binom_p = 0.5

# TUNING FOR THOMPSON ##################################################
m_thompson = 2

# MENZEL BOOTSTRAP ITERATIONS ##########################################
num_boot = 1000
# MNW BOOTSTRAP ITERATIONS #############################################
num_boot_MNW = 1000

# MAMEN RANDOM VARIABLE ################################################
mammen <- function(para){ a = para[1]; b = para[2]; p = para[3]; return( (p*a+(1-p)*b)^2 + (p*a^2+(1-p)*b^2-1)^2 + (p*a^3+(1-p)*b^3-1)^2 ); }
mammen_para = nlm(mammen,c(-1,1,0.5))$estimate
rmammen <- function(n,para){(para[1]-para[2])*(runif(n)<para[3])+para[2]}

# FUNCTION: EIGENVALUE-CORRECTION ######################################
EV_Correction <- function( A ){
 EV = eigen(A)
 Eval = EV$values
 Eval = (Eval >= 0) * Eval + (Eval < 0 ) * 0
 Evec = EV$vectors
 return( (Evec) %*% diag(Eval) %*% t(Evec) )
}

########################################################################
# BEGIN MONTE CARLO ####################################################
########################################################################
list_IID_95 = NULL
list_PNW_95 = NULL
list_CRi_95 = NULL
list_CRt_95 = NULL
list_CGM_95 = NULL
list_MNW_95 = NULL
list_MEN_95 = NULL
list_ME2_95 = NULL
list_THO_95 = NULL
list_DDG_95 = NULL
list_CHS_95 = NULL
for( iter in 1:NUM_MC ){

########################################################################
# DATA GENERATION (BEGIN)
########################################################################

### ALPHA (NORMAL) (N by 1) ############################################
alpha_0 = matrix(rnorm(N,0,1),N,1)
alpha_1 = matrix(rnorm(N,0,1),N,1)
alpha_2 = matrix(rnorm(N,0,1),N,1)
alpha_3 = matrix(rnorm(N,0,1),N,1)

### GAMMA (AR NORMAL) (1 by T) #########################################
if( DGP_AR ){
 gamma_0 = matrix(rnorm(1,0,1),1,T)
 gamma_1 = matrix(rnorm(1,0,1),1,T)
 gamma_2 = matrix(rnorm(1,0,1),1,T)
 gamma_3 = matrix(rnorm(1,0,1),1,T)
 for( idx in 2:T ){
  gamma_0[1,idx] = rho*gamma_0[1,idx-1] + rnorm(1,0,(1-rho^2)^0.5)
  gamma_1[1,idx] = rho*gamma_1[1,idx-1] + rnorm(1,0,(1-rho^2)^0.5)
  gamma_2[1,idx] = rho*gamma_2[1,idx-1] + rnorm(1,0,(1-rho^2)^0.5)
  gamma_3[1,idx] = rho*gamma_3[1,idx-1] + rnorm(1,0,(1-rho^2)^0.5)
 }
}

### GAMMA (MA NORMAL) (1 by T) #########################################
if( !DGP_AR ){
 pre_gamma_0 = matrix(rnorm(T+dim(MAcoef)[1]-1,0,1),,1)
 pre_gamma_1 = matrix(rnorm(T+dim(MAcoef)[1]-1,0,1),,1)
 pre_gamma_2 = matrix(rnorm(T+dim(MAcoef)[1]-1,0,1),,1)
 pre_gamma_3 = matrix(rnorm(T+dim(MAcoef)[1]-1,0,1),,1)
 gamma_0 = matrix(0,1,T)
 gamma_1 = matrix(0,1,T)
 gamma_2 = matrix(0,1,T)
 gamma_3 = matrix(0,1,T)
 for( idx in 1:T ){
  gamma_0[1,idx] = t(MAcoef) %*% matrix(pre_gamma_0[1:(dim(MAcoef)[1])+idx-1],,1)
  gamma_1[1,idx] = t(MAcoef) %*% matrix(pre_gamma_1[1:(dim(MAcoef)[1])+idx-1],,1)
  gamma_2[1,idx] = t(MAcoef) %*% matrix(pre_gamma_2[1:(dim(MAcoef)[1])+idx-1],,1)
  gamma_3[1,idx] = t(MAcoef) %*% matrix(pre_gamma_3[1:(dim(MAcoef)[1])+idx-1],,1)
 }
}

### EPSILON (NORMAL) (N by T) ##########################################
epsilon_0 = matrix(rnorm(N*T,0,1),N,T)
epsilon_1 = matrix(rnorm(N*T,0,1),N,T)

baseX = b[3,1] * ( alpha_1 %*% matrix(1,1,T) ) *
             ( matrix(1,N,1) %*% gamma_2 ) +
    b[4,1] * ( alpha_2 %*% matrix(1,1,T) ) *
             ( matrix(1,N,1) %*% gamma_1 ) +
    b[5,1] * epsilon_0

baseU = b[1,1] * ( alpha_0 %*% matrix(1,1,T) ) +
    b[2,1] * ( matrix(1,N,1) %*% gamma_0 ) +
    b[3,1] * ( alpha_1 %*% matrix(1,1,T) ) *
             ( matrix(1,N,1) %*% gamma_3 ) +
    b[4,1] * ( alpha_3 %*% matrix(1,1,T) ) *
             ( matrix(1,N,1) %*% gamma_1 ) +
    b[5,1] * epsilon_1

J = matrix(rbinom(N*T,binom_N,binom_p),N,T)

X = NULL
U = NULL
i_ = NULL
t_ = NULL
for( idx in 1:N ){
 for( tdx in 1:T ){
  X = c(X, baseX[idx,tdx] + b[6,1] * rnorm(J[idx,tdx],0,1))
  U = c(U, baseU[idx,tdx] + b[6,1] * rnorm(J[idx,tdx],0,1))
  i_ = c(i_ , array(idx,J[idx,tdx]))
  t_ = c(t_ , array(tdx,J[idx,tdx]))
 }
}

Y = true_beta[1,1] + true_beta[2,1] * X + U

### DOUBLE DIFFERENCE Y, X, AND U ######################################
meanYi = matrix(0,N,1)
meanXi = matrix(0,N,1)
meanUi = matrix(0,N,1)
for( idx in 1:N ){
 meanYi[idx] = mean(Y[i_==idx])
 meanXi[idx] = mean(X[i_==idx])
 meanUi[idx] = mean(U[i_==idx])
}
meanYt = matrix(0,T,1)
meanXt = matrix(0,T,1)
meanUt = matrix(0,T,1)
for( tdx in 1:T ){
 meanYt[tdx] = mean(Y[t_==tdx])
 meanXt[tdx] = mean(X[t_==tdx])
 meanUt[tdx] = mean(U[t_==tdx])
}
meanY = mean(Y)
meanX = mean(X)
meanU = mean(U)
for( idx in 1:N ){
 for( tdx in 1:T ){
  Y[i_==idx & t_==tdx] = Y[i_==idx & t_==tdx] - meanYi[idx] - meanYt[tdx] + meanY
  X[i_==idx & t_==tdx] = X[i_==idx & t_==tdx] - meanXi[idx] - meanXt[tdx] + meanX
  U[i_==idx & t_==tdx] = U[i_==idx & t_==tdx] - meanUi[idx] - meanUt[tdx] + meanU
}}

########################################################################
# DATA GENERATION (ENDED)
########################################################################



########################################################################
# ESTIMATION (BEGIN)
########################################################################

### OLS ################################################################
beta_hat = solve( t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) ) %*% t(cbind(1,matrix(c(X),,1))) %*% matrix(c(Y),,1)

### RESIDUALS ##########################################################
U_hat = Y - beta_hat[1,1] - beta_hat[2,1] * X

########################################################################
# ESTIMATION (ENDED)
########################################################################



########################################################################
# ESTIMATE IID OMEGA (BEGIN)
########################################################################
Omega_hat_iid = matrix(0,2,2)
Omega_hat_iid[1,1] = t(matrix(c(1*U_hat),,1)) %*% matrix(c(1*U_hat),,1)
Omega_hat_iid[2,1] = t(matrix(c(X*U_hat),,1)) %*% matrix(c(1*U_hat),,1)
Omega_hat_iid[1,2] = t(matrix(c(1*U_hat),,1)) %*% matrix(c(X*U_hat),,1)
Omega_hat_iid[2,2] = t(matrix(c(X*U_hat),,1)) %*% matrix(c(X*U_hat),,1)
Omega_hat_iid = 1/(N*T) * Omega_hat_iid
########################################################################
# ESTIMATE IID OMEGA (ENDED)
########################################################################

########################################################################
# ESTIMATE CHS OMEGA (BEGIN)
########################################################################

one_U_hat = matrix(0,N,T)
X_U_hat = matrix(0,N,T)
one_X = matrix(0,N,T)
for( idx in 1:N ){
 for( tdx in 1:T ){
  one_U_hat[idx,tdx] = sum( U_hat[i_==idx & t_==tdx] )
  X_U_hat[idx,tdx] = sum( X[i_==idx & t_==tdx] * U_hat[i_==idx & t_==tdx] )
  one_X[idx,tdx] = sum( X[i_==idx & t_==tdx] )
}}

### OMEGA HAT 1ST TERM #################################################
Omega_hat_1 = matrix(0,2,2)
for( i in 1:N ){
 Omega_hat_1[1,1] = Omega_hat_1[1,1] + sum( t(matrix(one_U_hat[i,],1,T)) %*% matrix(one_U_hat[i,],1,T) )
 Omega_hat_1[2,1] = Omega_hat_1[2,1] + sum( t(matrix(  X_U_hat[i,],1,T)) %*% matrix(one_U_hat[i,],1,T) )
 Omega_hat_1[1,2] = Omega_hat_1[1,2] + sum( t(matrix(one_U_hat[i,],1,T)) %*% matrix(  X_U_hat[i,],1,T) )
 Omega_hat_1[2,2] = Omega_hat_1[2,2] + sum( t(matrix(  X_U_hat[i,],1,T)) %*% matrix(  X_U_hat[i,],1,T) )
}
Omega_hat_1 = min(c(N,T))/N/(N*T^2) * Omega_hat_1

### OMEGA HAT 2ND TERM #################################################
Omega_hat_2 = matrix(0,2,2)
for( t in 1:T ){
 Omega_hat_2[1,1] = Omega_hat_2[1,1] + sum( matrix(one_U_hat[,t],N,1) %*% t(matrix(one_U_hat[,t],N,1)) )
 Omega_hat_2[2,1] = Omega_hat_2[2,1] + sum( matrix(  X_U_hat[,t],N,1) %*% t(matrix(one_U_hat[,t],N,1)) )
 Omega_hat_2[1,2] = Omega_hat_2[1,2] + sum( matrix(one_U_hat[,t],N,1) %*% t(matrix(  X_U_hat[,t],N,1)) )
 Omega_hat_2[2,2] = Omega_hat_2[2,2] + sum( matrix(  X_U_hat[,t],N,1) %*% t(matrix(  X_U_hat[,t],N,1)) )
}
Omega_hat_2 = min(c(N,T))/T/(N^2*T) * Omega_hat_2

### OMEGA HAT 3RD TERM #################################################
Omega_hat_3 = matrix(0,2,2)
Omega_hat_3[1,1] = t(matrix(c(one_U_hat),,1)) %*% matrix(1*c(one_U_hat),,1)
Omega_hat_3[2,1] = t(matrix(c(  X_U_hat),,1)) %*% matrix(1*c(one_U_hat),,1)
Omega_hat_3[1,2] = t(matrix(c(one_U_hat),,1)) %*% matrix(1*c(  X_U_hat),,1)
Omega_hat_3[2,2] = t(matrix(c(  X_U_hat),,1)) %*% matrix(1*c(  X_U_hat),,1)
Omega_hat_3 = min(c(N,T))/(N*T)/(N*T) * Omega_hat_3

### OMEGA HAT 4TH TERM FOR THOMPSON ####################################
Omega_hat_4_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_4_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  temp_Omega_hat_4_thompson[1,1] = temp_Omega_hat_4_thompson[1,1] + sum( matrix(one_U_hat[,t],N,1) %*% t(matrix(one_U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[2,1] = temp_Omega_hat_4_thompson[2,1] + sum( matrix(  X_U_hat[,t],N,1) %*% t(matrix(one_U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[1,2] = temp_Omega_hat_4_thompson[1,2] + sum( matrix(one_U_hat[,t],N,1) %*% t(matrix(  X_U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[2,2] = temp_Omega_hat_4_thompson[2,2] + sum( matrix(  X_U_hat[,t],N,1) %*% t(matrix(  X_U_hat[,t-j],N,1)) )
 }
 temp_Omega_hat_4_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_4_thompson
 Omega_hat_4_thompson = Omega_hat_4_thompson + temp_Omega_hat_4_thompson
}
Omega_hat_4_thompson = min(c(N,T))/T * Omega_hat_4_thompson

### OMEGA HAT 5TH TERM FOR THOMPSON ####################################
Omega_hat_5_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_5_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  temp_Omega_hat_5_thompson[1,1] = temp_Omega_hat_5_thompson[1,1] + sum( matrix(one_U_hat[,t-j],N,1) %*% t(matrix(one_U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[2,1] = temp_Omega_hat_5_thompson[2,1] + sum( matrix(  X_U_hat[,t-j],N,1) %*% t(matrix(one_U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[1,2] = temp_Omega_hat_5_thompson[1,2] + sum( matrix(one_U_hat[,t-j],N,1) %*% t(matrix(  X_U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[2,2] = temp_Omega_hat_5_thompson[2,2] + sum( matrix(  X_U_hat[,t-j],N,1) %*% t(matrix(  X_U_hat[,t],N,1)) )
 }
 temp_Omega_hat_5_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_5_thompson
 Omega_hat_5_thompson = Omega_hat_5_thompson + temp_Omega_hat_5_thompson
}
Omega_hat_5_thompson = min(c(N,T))/T * Omega_hat_5_thompson

### OMEGA HAT 6TH TERM FOR THOMPSON ####################################
Omega_hat_6_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_6_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  temp_Omega_hat_6_thompson[1,1] = temp_Omega_hat_6_thompson[1,1] + sum( matrix(one_U_hat[,t],N,1) * (matrix(one_U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[2,1] = temp_Omega_hat_6_thompson[2,1] + sum( matrix(  X_U_hat[,t],N,1) * (matrix(one_U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[1,2] = temp_Omega_hat_6_thompson[1,2] + sum( matrix(one_U_hat[,t],N,1) * (matrix(  X_U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[2,2] = temp_Omega_hat_6_thompson[2,2] + sum( matrix(  X_U_hat[,t],N,1) * (matrix(  X_U_hat[,t-j],N,1)) )
 }
 temp_Omega_hat_6_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_6_thompson
 Omega_hat_6_thompson = Omega_hat_6_thompson + temp_Omega_hat_6_thompson
}
Omega_hat_4_thompson = min(c(N,T))/T * Omega_hat_4_thompson

### OMEGA HAT 7TH TERM FOR THOMPSON ####################################
Omega_hat_7_thompson = matrix(0,2,2)
for( j in 1:m_thompson ){
 temp_Omega_hat_7_thompson = matrix(0,2,2)
 for( t in (j+1):T ){
  temp_Omega_hat_7_thompson[1,1] = temp_Omega_hat_7_thompson[1,1] + sum( matrix(one_U_hat[,t-j],N,1) * (matrix(one_U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[2,1] = temp_Omega_hat_7_thompson[2,1] + sum( matrix(  X_U_hat[,t-j],N,1) * (matrix(one_U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[1,2] = temp_Omega_hat_7_thompson[1,2] + sum( matrix(one_U_hat[,t-j],N,1) * (matrix(  X_U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[2,2] = temp_Omega_hat_7_thompson[2,2] + sum( matrix(  X_U_hat[,t-j],N,1) * (matrix(  X_U_hat[,t],N,1)) )
 }
 temp_Omega_hat_7_thompson = (j<=m_thompson) / (N^2*T) * temp_Omega_hat_7_thompson
 Omega_hat_7_thompson = Omega_hat_7_thompson + temp_Omega_hat_7_thompson
}
Omega_hat_7_thompson = min(c(N,T))/T * Omega_hat_7_thompson

### COMPUTE M HAT ######################################################
S = colSums(X_U_hat)[2:T]
Slag = colSums(X_U_hat)[1:(T-1)]
rho_hat = cov(S,Slag)/var(Slag)
m = 1.8171 * (rho_hat^2/(1-rho_hat^2)^2)^(1/3) * T^(1/3)
m = (min(c(m,T-1)))

### OMEGA HAT 4TH TERM WITH BARTLET KERNEL #############################
Omega_hat_4 = matrix(0,2,2)
if( m >= 1 ){
 for( j in 1:m ){
  temp_Omega_hat_4 = matrix(0,2,2)
  w = 1 - (j/(m+1))
  for( t in (j+1):T ){
   temp_Omega_hat_4[1,1] = temp_Omega_hat_4[1,1] + sum( matrix(one_U_hat[,t],N,1) %*% t(matrix(one_U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[2,1] = temp_Omega_hat_4[2,1] + sum( matrix(  X_U_hat[,t],N,1) %*% t(matrix(one_U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[1,2] = temp_Omega_hat_4[1,2] + sum( matrix(one_U_hat[,t],N,1) %*% t(matrix(  X_U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[2,2] = temp_Omega_hat_4[2,2] + sum( matrix(  X_U_hat[,t],N,1) %*% t(matrix(  X_U_hat[,t-j],N,1)) )
  }
  temp_Omega_hat_4 = w / (N^2*T) * temp_Omega_hat_4
  Omega_hat_4 = Omega_hat_4 + temp_Omega_hat_4
 }
 Omega_hat_4 = min(c(N,T))/T * Omega_hat_4
}

### OMEGA HAT 5TH TERM WITH BARTLET KERNEL #############################
Omega_hat_5 = matrix(0,2,2)
if( m >= 1 ){
 for( j in 1:m ){
  temp_Omega_hat_5 = matrix(0,2,2)
  w = 1 - (j/(m+1))
  for( t in (j+1):T ){
   temp_Omega_hat_5[1,1] = temp_Omega_hat_5[1,1] + sum( matrix(one_U_hat[,t-j],N,1) %*% t(matrix(one_U_hat[,t],N,1)) )
   temp_Omega_hat_5[2,1] = temp_Omega_hat_5[2,1] + sum( matrix(  X_U_hat[,t-j],N,1) %*% t(matrix(one_U_hat[,t],N,1)) )
   temp_Omega_hat_5[1,2] = temp_Omega_hat_5[1,2] + sum( matrix(one_U_hat[,t-j],N,1) %*% t(matrix(  X_U_hat[,t],N,1)) )
   temp_Omega_hat_5[2,2] = temp_Omega_hat_5[2,2] + sum( matrix(  X_U_hat[,t-j],N,1) %*% t(matrix(  X_U_hat[,t],N,1)) )
  }
  temp_Omega_hat_5 = w / (N^2*T) * temp_Omega_hat_5
  Omega_hat_5 = Omega_hat_5 + temp_Omega_hat_5
 }
 Omega_hat_5 = min(c(N,T))/T * Omega_hat_5
}

########################################################################
# ESTIMATE CHS OMEGA (BEGIN)
########################################################################



########################################################################
# VARIANCE ESTIMATION (BEGIN)
########################################################################

### DESIGN MATRIX Q HAT ################################################
Q_hat = t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) / (N*T)

### STANDARD ERRORS ####################################################
SE_IID = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_iid )                                                                                                                       )%*%solve(Q_hat)/(N*T) ) )
SE_CRi = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 )                                                                                                                         )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CRt = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_2 )                                                                                                                         )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CGM = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 )                                                                                             )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_DDG = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 )                                                                                                           )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_THO = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 + Omega_hat_4_thompson + Omega_hat_5_thompson - Omega_hat_6_thompson - Omega_hat_7_thompson ) )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CHS = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 + Omega_hat_4 + Omega_hat_5 )                                                                 )%*%solve(Q_hat)/min(c(N,T)) ) )

########################################################################
# VARIANCE ESTIMATION (ENDED)
########################################################################



########################################################################
# WILD BOOTSTRAP CRITICAL VALUE (MNW2021) (BEGIN)
########################################################################
#MNW_t_list = NULL
#for( boot_iterations in 1:num_boot_MNW ){
#rademacher = -1+2*(runif(length(U_hat))<0.5)
#U_hat_wild = rademacher * U_hat
#Y_boot = beta_hat[1,1] + beta_hat[2,1] * X + U_hat_wild
#beta_hat_boot = solve( t(cbind(1,matrix(c(X),,1))) %*% cbind(1,matrix(c(X),,1)) ) %*% t(cbind(1,matrix(c(X),,1))) %*% matrix(c(Y_boot),,1)
#U_hat_boot = Y_boot - beta_hat_boot[1,1] - beta_hat_boot[2,1] * X
#
#one_U_hat_boot = matrix(0,N,T)
#X_U_hat_boot = matrix(0,N,T)
#one_X_boot = matrix(0,N,T)
#for( idx in 1:N ){
# for( tdx in 1:T ){
#  one_U_hat_boot[idx,tdx] = sum( U_hat_boot[i_==idx & t_==tdx] )
#  X_U_hat_boot[idx,tdx] = sum( X[i_==idx & t_==tdx] * U_hat_boot[i_==idx & t_==tdx] )
#  one_X_boot[idx,tdx] = sum( X[i_==idx & t_==tdx] )
#}}
#
### OMEGA HAT 1ST TERM #################################################
#Omega_hat_1_boot = matrix(0,2,2)
#for( i in 1:N ){
# Omega_hat_1_boot[1,1] = Omega_hat_1_boot[1,1] + sum( t(matrix(one_U_hat_boot[i,],1,T)) %*% matrix(one_U_hat_boot[i,],1,T) )
# Omega_hat_1_boot[2,1] = Omega_hat_1_boot[2,1] + sum( t(matrix(  X_U_hat_boot[i,],1,T)) %*% matrix(one_U_hat_boot[i,],1,T) )
# Omega_hat_1_boot[1,2] = Omega_hat_1_boot[1,2] + sum( t(matrix(one_U_hat_boot[i,],1,T)) %*% matrix(  X_U_hat_boot[i,],1,T) )
# Omega_hat_1_boot[2,2] = Omega_hat_1_boot[2,2] + sum( t(matrix(  X_U_hat_boot[i,],1,T)) %*% matrix(  X_U_hat_boot[i,],1,T) )
#}
#Omega_hat_1_boot = min(c(N,T))/N/(N*T^2) * Omega_hat_1_boot
#
### OMEGA HAT 2ND TERM #################################################
#Omega_hat_2_boot = matrix(0,2,2)
#for( t in 1:T ){
# Omega_hat_2_boot[1,1] = Omega_hat_2_boot[1,1] + sum( matrix(one_U_hat_boot[,t],N,1) %*% t(matrix(one_U_hat_boot[,t],N,1)) )
# Omega_hat_2_boot[2,1] = Omega_hat_2_boot[2,1] + sum( matrix(  X_U_hat_boot[,t],N,1) %*% t(matrix(one_U_hat_boot[,t],N,1)) )
# Omega_hat_2_boot[1,2] = Omega_hat_2_boot[1,2] + sum( matrix(one_U_hat_boot[,t],N,1) %*% t(matrix(  X_U_hat_boot[,t],N,1)) )
# Omega_hat_2_boot[2,2] = Omega_hat_2_boot[2,2] + sum( matrix(  X_U_hat_boot[,t],N,1) %*% t(matrix(  X_U_hat_boot[,t],N,1)) )
#}
#Omega_hat_2_boot = min(c(N,T))/T/(N^2*T) * Omega_hat_2_boot
#
### OMEGA HAT 3RD TERM #################################################
#Omega_hat_3_boot = matrix(0,2,2)
#Omega_hat_3_boot[1,1] = t(matrix(c(one_U_hat_boot),,1)) %*% matrix(1*c(one_U_hat_boot),,1)
#Omega_hat_3_boot[2,1] = t(matrix(c(  X_U_hat_boot),,1)) %*% matrix(1*c(one_U_hat_boot),,1)
#Omega_hat_3_boot[1,2] = t(matrix(c(one_U_hat_boot),,1)) %*% matrix(1*c(  X_U_hat_boot),,1)
#Omega_hat_3_boot[2,2] = t(matrix(c(  X_U_hat_boot),,1)) %*% matrix(1*c(  X_U_hat_boot),,1)
#Omega_hat_3_boot = min(c(N,T))/(N*T)/(N*T) * Omega_hat_3_boot
#
#SE_CGM_boot = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1_boot + Omega_hat_2_boot - Omega_hat_3_boot ) )%*%solve(Q_hat)/min(c(N,T)) ) )
#t_wild_boot = (beta_hat_boot - beta_hat ) / matrix(SE_CGM_boot,2,1)
#MNW_t_list = cbind( MNW_t_list, t_wild_boot )
#}
#crit_MNW = apply(abs(MNW_t_list),1,quantile,p=0.95)
########################################################################
# WILD BOOTSTRAP CRITICAL VALUE (MNW2021) (ENDED)
########################################################################



cov_IID_95 = t(abs(beta_hat - true_beta) <= SE_IID*qnorm(0.975)); if( is.na(cov_IID_95)[1] ){cov_IID_95[1]=0}; if( is.na(cov_IID_95)[2] ){cov_IID_95[2]=0};
cov_CRi_95 = t(abs(beta_hat - true_beta) <= SE_CRi*qnorm(0.975)); if( is.na(cov_CRi_95)[1] ){cov_CRi_95[1]=0}; if( is.na(cov_CRi_95)[2] ){cov_CRi_95[2]=0};
cov_CRt_95 = t(abs(beta_hat - true_beta) <= SE_CRt*qnorm(0.975)); if( is.na(cov_CRt_95)[1] ){cov_CRt_95[1]=0}; if( is.na(cov_CRt_95)[2] ){cov_CRt_95[2]=0};
cov_CGM_95 = t(abs(beta_hat - true_beta) <= SE_CGM*qnorm(0.975)); if( is.na(cov_CGM_95)[1] ){cov_CGM_95[1]=0}; if( is.na(cov_CGM_95)[2] ){cov_CGM_95[2]=0};
#cov_MNW_95 = t(abs(beta_hat - true_beta) <= SE_CGM*crit_MNW);     if( is.na(cov_CGM_95)[1] ){cov_CGM_95[1]=0}; if( is.na(cov_CGM_95)[2] ){cov_CGM_95[2]=0};
cov_THO_95 = t(abs(beta_hat - true_beta) <= SE_THO*qnorm(0.975)); if( is.na(cov_THO_95)[1] ){cov_THO_95[1]=0}; if( is.na(cov_THO_95)[2] ){cov_THO_95[2]=0};
cov_DDG_95 = t(abs(beta_hat - true_beta) <= SE_DDG*qnorm(0.975)); if( is.na(cov_DDG_95)[1] ){cov_DDG_95[1]=0}; if( is.na(cov_DDG_95)[2] ){cov_DDG_95[2]=0};
cov_CHS_95 = t(abs(beta_hat - true_beta) <= SE_CHS*qnorm(0.975)); if( is.na(cov_CHS_95)[1] ){cov_CHS_95[1]=0}; if( is.na(cov_CHS_95)[2] ){cov_CHS_95[2]=0};

list_IID_95 = rbind(list_IID_95, cov_IID_95)
list_CRi_95 = rbind(list_CRi_95, cov_CRi_95)
list_CRt_95 = rbind(list_CRt_95, cov_CRt_95)
list_CGM_95 = rbind(list_CGM_95, cov_CGM_95)
#list_MNW_95 = rbind(list_MNW_95, cov_MNW_95)
list_THO_95 = rbind(list_THO_95, cov_THO_95)
list_DDG_95 = rbind(list_DDG_95, cov_DDG_95)
list_CHS_95 = rbind(list_CHS_95, cov_CHS_95)

options(digits=3)
print(c(iter,colMeans(list_IID_95)[2],colMeans(list_CRi_95)[2],colMeans(list_CRt_95)[2],colMeans(list_CGM_95)[2],colMeans(list_THO_95)[2],colMeans(list_CHS_95)[2]))
Sys.sleep(0.00000000000000001); flush.console()
}
########################################################################
# END OF MONTE CARLO ###################################################
########################################################################


RESULTS = cbind(colMeans(list_IID_95),colMeans(list_CRi_95),colMeans(list_CRt_95),colMeans(list_CGM_95),colMeans(list_THO_95),colMeans(list_CHS_95))
colnames(RESULTS) <- c("EHW_95","CRi_95","CRt_95","CGM_95","Tho_95","CHS_95")
rownames(RESULTS) <- c("beta0","beta1")
RESULTS
if(DGP_AR){
 filename = paste("results_AR_f_",b[1],"_",b[2],"_",b[3],"_rho",rho,"_N",N,"_T",T,".csv",sep="")
}else{
 filename = paste("results_MA_f_",b[1],"_",b[2],"_",b[3],"_q",q,"_N",N,"_T",T,".csv",sep="")
}
write.csv(RESULTS,filename)
